Unjamming of Granular Packings due to Local Perturbations: Stability and Decay of 

Displacements 
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We study the mechanical response generated by local deformations in jammed packings of rigid 
disks. Based on discrete element simulations we determine the critical force of the local perturbation 
that is needed to break the mechanical equilibrium and examine the generated displacement field. 
Displacements decay as a power law of the distance from the perturbation point. The decay exponent 
and the critical force exhibit nontrivial dependence on the friction: Both quantities are nonmonotonic 
and have a sharp maximum at the friction coefficient 0.1. We find that the mechanical response 
properties are closely related to the problem of force-indeterminacy where similar nonmonotonic 
behavior was observed previously. We establish direct connection between the critical force and the 
ensemble of static force networks. 
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Granular materials constitute an ideal field for study- 
ing the physics of the jamming transition [l], 0]. One 
of the most exciting challenges of research is to provide 
better understanding of the onset of yielding in granu- 
lar media, an example of unjamming. When the external 
load on a static assembly of grains is changed at a certain 
point the load may become incompatible with the inner 
structure of the packing and the solid state looses its sta- 
bility. How exactly this happens on grain-scale and what 
are the key features of the transition between statics and 
flow are intriguing and unresolved problems [H, 0, [H, @ . 
Moreover, it is of essential importance in many appli- 
cations to be able to predict, initiate or prevent such 
transitions. 

In this Letter we study yielding induced by local per- 
turbations and focus on two response properties. First, 
we investigate the resistance of the packing that is ex- 
erted against a local deformation. Second, based on the 
generated displacement field, we study how far the ef- 
fect of the perturbation penetrates into the packing. In 
recent years local perturbations were extensively used to 
study the nature of stress transmission in granular media 
@> H Q ■ These perturbations were typically weak in the 
sense that they did not break the stability: The origi- 
nal contact network can maintain the perturbation, par- 
ticle displacements are reversible and occur merely due 
to elastic distortions of contacts [8[. Here we deal with 
rigid (undcformable) grains with help of the simulation 
technique of contact dynamics [1(3, El therefore elastic 
deformations are excluded. As opposed to the weak per- 
turbations, in our case, plastic deformation of the packing 
is initiated 0, Q ■ Plastic rearrangements caused by local 
perturbations have been studied in recent experiments [5[ 
where interesting power law decay of the rearrangement 
field was found. 

An important question we address here is what effect 
the particle-particle friction fi has on the response of the 



packings. This was motivated by the prediction that re- 
sponse properties may show nonmonotonic behavior as 
the function of /i, i.e. small and large frictions lead to 
similar behavior which is different from the behavior for 
intermediate friction [T^]. This idea was based on the 
ensemble of admissible force configurations. 

It is known that, in general case, mechanical equilib- 
rium and Coulomb condition do not determine contact 
forces uniquely 0, 0, [T3, EH • Consequently there is an 
ensemble of force networks that satisfy these conditions 
in the same contact geometry and for the same exter- 
nal load. It depends strongly on friction how large the 
indeterminacy of individual contact forces is. E.g. for 
frictionless rigid grains the indeterminacy vanishes and 
the ensemble shrinks to a single admissible force network 
[T3 |. It was found for 2D packings of disks that inde- 
terminacy of forces becomes the largest for [i ss 0.1 
further increase of the friction leads again to smaller in- 
determinacy. It was argued in [12| that the extent of 
force-indeterminacy must affect the stability of the pack- 
ing, because the more freedom the forces have the more 
possibility the packing has to resolve changes of the load 
without rearrangements. This suggests that mechanical 
response properties may show nonmonotonic dependence 
on the coefficient of friction. The main message of the 
present work is that such correlation indeed exists. 

The connection between mechanical properties and 
force indeterminacy was also examined from an other 
point of view in [14j . where the maximum admissible 
shear tress was deduced from the ensemble of force net- 
works for packings of frictionless soft particles. 

The systems we examine are two-dimensional random 
packings of perfectly rigid disks in zero gravity. The unit 
of the length is set to the maximum grain radius, radii 
are distributed uniformly between 0.5 and 1. The num- 
ber of the grains contained by the packings ranges from 
500 to 8000. Our numerical experiments consist of two 
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parts. First by compression we generate dense random 
packings then we probe the packings by perturbing sin- 
gle contacts. We apply the contact dynamics method for 
both procedures. 

The compaction starts from a gas-like state with grains 
randomly distributed in a square-shaped cell. Periodic 
boundary conditions are applied in both directions and 
there are no walls in the system. Instead of using 
pistons the compaction is achieved by the method of 
Andersen [lH]: We impose constant external pressure p ex t 
and let the volume of the cell evolve in time. A detailed 
description of the coupling between the pressure bath and 
the system can be found in [l5j . Here we mention only 
the favorable feature of this method that it gives rise to 
homogeneous compression and, due to exclusion of side 
walls, boundary effects are avoided. 

As the size of the cell decreases, grains form contacts 
and start building up the inner pressure in order to avoid 
interpenctration. Finally a static, jammed configuration 
is reached where the grains block further compaction. 
In the final packing contact forces are such that they 
provide mechanical equilibrium for each grain and the 
corresponding inner pressure p- ln equals to p eK t ■ 

After that we turn to the perturbation part where we 
choose a pair of contacting grains and force them to move 
apart. One way of doing this would be to apply a small 
normal force between the two grains and continually in- 
crease it. It is expected that if the perturbation force 
is small enough then it can be resolved by the packing 
without rearrangements. If the force is increased further 
the yield point will be reached where the perturbation 
induce sliding and/or opening of some contacts and ini- 
tiates collective rearrangements of the particles at least 
in the vicinity of the perturbation point. 

It is very time consuming in computer simulations to 
scan a region of the perturbation force in order to find 
the yield point. Therefore we apply another method. In- 
stead of tuning the force we control the deformation. The 
idea is to bring the system immediately to the yield point 
by enforcing the opening of the contact: We prescribe a 
small gap at the contact of perturbation and then de- 
termine the force that is needed to fulfill this constraint. 
This concept is suited very well to the contact dynam- 
ics method where interparticle forces are handled as con- 
straint forces, i.e. they are calculated based on constraint 
conditions which prescribe the relative motion of the con- 
tact surfaces |Tl| . It is beneficial to choose small gap 
size as we are interested in the onset of motion, how the 
static structure breaks due to the perturbation. Large de- 
formations, e.g. creation of new contacts, are out of the 
scope of the present study. We checked that for small 
gap sizes the displacement field (up to a constant factor) 
and the critical perturbation force become independent 
of the size of the gap. Our numerical measurements are 
performed in this region, the size of the gap g is set to 
10- 9 . 

With the above technique we open up the contact that 
is selected for the perturbation and measure the gener- 
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FIG. 1: Displacement response field in the contact frame av- 
eraged over several thousand perturbations. The system con- 
tained 3000 disks with friction coefficient 0.5. 

ated displacement field and the required perturbation 
force. The latter describes the strength of the system 
against the local perturbation and we will refer to it as 
the critical force F cr ; t . Both the force and the displace- 
ment response depend strongly on the place of the per- 
turbation. First we check their average properties. 

By studying the displacement field our goal is to find 
out how far the rearrangements have to penetrate into 
the packing in order to allow the prescribed local defor- 
mation. Is there a related length scale? When a single 
contact is perturbed the displacements of particle centers 
form a disordered vector field which can be widespread 
or more localized depending on the perturbed contact. 
In order to calculate the average displacement field we 
perturb all contacts one by one always starting with the 
original static packing. In each case the particle move- 
ments are recorded in the local contact frame where the 
perturbed contact sits in the origin and the x-axis is cho- 
sen parallel to the contact normal, i.e. x indicates the di- 
rection of the separation. Fig. [T] shows a smooth displace- 
ment field obtained by averaging over the perturbed con- 
tacts. The apparent circular shape of the system is the 
consequence of the averaging because the original square 
shape can have different orientations when transformed 
into different contact frames. 

The quadrupolar structure in Fig. [T] is interestingly 
very similar to the displacement fields that have been 
observed in sheared systems of deformable frictionless 
grains (l6j , where localized quadrupolar deformations ap- 
pear at the onset of plastic events. Here, for further anal- 
ysis we consider only the magnitude of the displacement 
vectors and their distance r from the place of the pertur- 
bation (angle of the position is averaged out). The decay 
of the the average displacement d is shown in Fig. [2j 

The curves d(r) show power law decay 

d(xr- a . (1) 
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FIG. 2: The magnitude of the displacements, d, in terms 
of the distance from the perturbed contact, r. g stands for 
the gap generated at the perturbed contact. Different slopes 
correspond to different friction coefficients \i. For each value 
of friction four systems of different sizes are investigated. The 
total number of particles are between 500 and 8000. 



Thus we find that no penetration length can be related 
to the rearrangements, but instead, the penetration can 
be best characterized by an exponent (a). 

Variation of the system size (Fig. [2]), within the modest 
range allowed by the numerical tools, shows no system- 
atic change in the exponent. There is, however, a strong 
dependence on friction. The values we obtained for a 
are between 0.7 and 1.4 (T7| . Similar power law behavior 
with the same range of a has been found experimentally 
by Kolb et al.[| by moving an intruder in a system of 
disks. Next we discuss the influence of friction on the 
exponent and on the critical force. 

It is important to note that the packings, that are 
subjected to perturbations, are newly generated for each 
value of friction therefore different values of /i are ac- 
companied by different packing structures. In Fig. [3ja 
the average contact number z [18| and the average con- 
tact force (fcont} is plotted. In F cont we take only the 
normal component of the contact force into account, the 
average ( ) is meant over all the contacts of the given 
packing. The contact forces are measured in units F 
set by the external pressure and by the average radius of 
the particles, Fq = 2R avg P cxt - The behavior of the pack- 
ing fraction (not shown) is similar to that of z, its value 
is 0.84 (0.80) in the low (high) friction limit. Quantities, 
that describe the properties of the packing, have the com- 
mon behavior that they exhibit plateaus for low and high 
friction and show a smooth and monotonic transition in 
between. 

We find a completely different behavior in the response 
(Fig. Ob). Both a and the average critical force (F cr i t ) 
are nonmonotonic functions of the friction coefficient and 
exhibit a quite sharp peak at the same place /x ~ 0.1. 
This is a remarkable behavior because all systems are 
subjected to exactly the same compaction and perturba- 
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FIG. 3: (a) Influence of friction on the average normal contact 
force (F cont ) (full circles) and on the contact number z (open 
circles), (b) Average critical force (F cr it) (full circles) with re- 
spect to the average normal contact force and the penetration 
exponent a (open circles) as functions of the friction coeffi- 
cient fi. The horizontal lines show the reference value I for 
the exponent a (dashed) and for the critical force (solid). The 
inset shows the role of the structural change in the packing. 
For normally constructed packings (F cr it) is plotted with full 
squares (new packing is constructed for each friction). Open 
squares were obtained by changing the friction in a fixed pack- 
ing configuration. 



tion procedures, friction was the only quantity that has 
been changed. 

When the friction is increased starting from zero, pack- 
ings are getting stronger against the perturbation and 
the induced rearrangements become more localized. At 
fj, ps 0.1 the process takes a sharp turn and further in- 
crease of the friction leads to softening and derealization. 

The nonmonotonic behavior and the position of the 
maxima in (Fig. [3]b) are reminiscent of the force- 
indeterminacy that was reported in [r3 |. In order to 
clarify what role the force- indeterminacy plays in the me- 
chanical response we perform two independent numerical 
measurements in the same test system: on the one hand 
we investigate the force-indeterminacy on the other hand 
we determine the critical forces. Then the data are com- 
pared on the level of single contacts. 

As mentioned before, in frictional case contact forces 
are not determined uniquely but there is an ensemble 
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FIG. 4: The critical force F cr it (cross), the original contact 
force (open circle) and the possible values of the normal con- 
tact forces (dots, mostly merged to intervals) are shown for 
each contact of the packing. 



of force networks that correspond to equilibrium and to 
the compatibility conditions (forces are limited by the 
Coulomb cone). This ensemble of force networks forms 
a convex set in the space of the contact forces [l^ . We 
explore this high dimensional solution set using a ran- 
dom walk starting with the original force network. The 
exploration process is the same as in [li^ . Each step 
of the random walk provides one possible force-network. 
Because the dimension of the solution set is large (pro- 
portional to the number of contacts) only small systems 
can be handled efficiently. We generate a packing of 50 
disks with \i = 0.1 (total number of contacts is 85). As a 
consequence of the convexity, the possible values of the 
normal (the tangential) component of any contact force 
are defined by an interval. These intervals are traced 
out in Fig. fj] with help of 20 000 force networks that are 
collected by the exploration procedure. 

The figure shows also the original contact forces (pro- 



vided by the compaction process) and the critical forces. 
The data reveal that the critical force obtained by the 
perturbation coincides with the maximum of the equilib- 
rium solutions at each contact. Thus the critical force 
can be directly obtained from the force-ensemble. It also 
can be seen that a strong (weak) contact has typically 
large (small) critical force. 

This picture also explains the limit of small and large 
frictions where the indeterminacy vanishes. In these 
cases the length of the force-intervals goes to zero. Which 
means that a pair of contacting particles can not resist a 
force of separation larger than the force itself that orig- 
inally presses the two contact surfaces together. This is 
why the average critical force approaches (-Fcont) on the 
left and right side of Fig. Ob. (F crit ) can differ signifi- 
cantly from (-Fcont) only if the indeterminacy of forces is 
large in the packing. 

It was argued in [12j that the nonmonotonic force- 
indeterminacy can be attributed to two competing effect 
of the increasing friction: first, the strengthening of the 
contacts (increase), and second, the structural change 
inside the packing (decrease) . Accordingly, if the second 
effect is switched off the decay part of the curve -F C rit(/-t) 
disappears as shown by the inset of Fig. [3Jb. Here the 
response of the same packing configuration was tested for 
various values of friction. 

In this Letter we reported nontrivial dependence of the 
mechanical response on the coefficient of friction and re- 
lated the results to the indeterminacy of forces. Nonethe- 
less, many questions arise that urge further studies. We 
list three of them here. What is the origin of the dis- 
placement exponent a? Numerical data show strong cor- 
relation between a and the critical force. What is the 
connection between them? Finally, how does the pic- 
ture, which is presented here for rigid particles, change 
for deformable ones? 
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